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1.  INTRODUCTION 


1.1  THE  FINITE  ELEMENT  METHOD 

The  finite  element  method,  in  general,  is  an  approximate  method  to 
solve  differential  equations.  Using  variational  calculus  the  differential 
equation  under  consideration  is  posed  as  a  functional.  The  resulting 
functional  depends  upon  the  unknowns  and  their  derivatives  with  respect  to 
the  spatial  coordinates  x,y  and  z  and  possibly  the  time,  t.  In  structural 
problems  the  functional  represents  a  meaningful  quantity,  namely,  the 
potential  energy.  However,  in  general,  the  functional  may  not  have  any 
physical  interpretation.  Minimizing  the  functional  with  respect  to  the 
unknowns  is  equivalent  to  solving  the  differential  equation.  The  functional 
is  minimized  by  setting  its  first  variation  to  zero.  In  structural  problems  this 
corresponds  to  the  well  known  concept  of  minimization  of  the  potential 
energy.  The  result  of  the  minimization  is  a  set  of  algebraic  equations 

(1.1)  [  K  ]{  u  }  =  {f} 

where  [  K  ]  is  the  matrix  of  coefficients  of  the  unknowns  and  is  known  as  the 

*  „  l  r  r  .  ,  ,  i  ,  r .  „ 
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"stiffness  matrix", 

{  u  }  is  an  array  of  unknowns  and 
{  f  }  is  an  array  of  forcing  functions. 

The  equations  are  then  solved  for  u.  In  a  typical  application  the 
domain  under  consideration  is  modeled  by  dividing  it  into  elements.  An 
interpolation  function,  or  shape  function,  is  set  for  the  elements  to  interpolate 
values  of  the  unknown  at  any  point  inside  the  element  in  terms  of  its  values 
at  the  nodes.  This  interpolation  function  is  used  in  the  functional  which  when 
minimized  as  described  above,  yields  the  stiffness  matrix.  The  reader  is 
directed  to  several  excellent  books  on  finite  element  method  by  Zienkiewics 
[1,2]*,  Segerlind  [3],  Reddy  [4],  and  Huston  [5].  The  point  to  note  for  this 
report  is  the  important  role  of  the  minimization  process  involved  in  the  finite 
element  methods. 


1 2  MESH  GENERATION 

Each  element  in  the  finite  element  model  is  addressed  by  its 


•Numbers  in  brackets  indicate  references  listed  at  the  end. 
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number.  Also  each  node  is  addressed  by  its  number.  The  inter-connectivity 
of  the  elements  is  determined  by  the  common  nodes  shared  by  the 
elements.  In  a  model  with  few  elements  and  nodes,  the  user  can  manually 
divide  the  domain,  number  each  element  and  node,  and  keep  track  of  the 
element  connectivity.  However,  in  models  with  many  nodes  and  elements, 
the  effort  required  to  divide  the  domain  into  elements  and  attend  to 
connectivity  is  great.  It  then  becomes  difficult  to  accomplish  this  task 
without  committing  errors.  However,  there  are  several  finite  element  pre¬ 
processors  which  do  this  job  automatically  once  the  geometry  is  defined. 
Users  can  then  devote  more  time  to  interpreting  results.  Shephard  [6] 
has  reviewed  the  current  trends  in  mesh  generation.  Although  there  are 
several  ways  to  generate  meshes,  these  methods  fall  into  two  broad 
categories: 

1.2.1  MAPPING  TECHNIQUES 

This  type  of  mesh  generation  is  best  suited  when  the  geometry  is 
simple  -  as  in  the  case  of  a  rectangle  or  a  cuboid.  Typically  the  user 
needs  to  choose  the  number  of  elements  on  each  of  the  edges  that 
defines  the  geometry  and  the  element  concentration  along  the  edges.  The 
software  then  generates  the  mesh  simply  by  joining  nodes  on  the  opposite 
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edges.  NASTRAN  MSGMESHt  ,  GIFTSn  ,  SUPERSAP*  and 
SUPARTABf  (in  I-DEAS)  have  this  capability.  For  a  more  complicated 
geometry  Schwarz-Christoffel  [7]  mapping  has  been  used  .  The  difficulty  in 
evaluating  integrals  involved  in  the  Schwarz-Christoffel  transformation 
however  makes  this  technique  less  attractive.  Moreover,  mesh  generated  by 
these  techniques  may  introduce  elements  with  high  aspect  ratios  and 
elements  that  are  highly  distorted. 


1.2.2  FREE  MESH  GENERATION 

This  method  of  generation  is  best  suited  for  models  with 
complicated  geometry.  SUPERTAB  has  this  capability.  The  model  is 
broken  down  into  sub-areas  and  sub-volumes.  On  each  of  the  curves  of 
every  sub-area  and  sub-volume  the  number  of  elements  and  their 
concentrations  are  selected.  The  software  then  generates  a  mesh  that  is 
consistent  with  the  selected  values  and  satisfies  the  requirement  on  the 
aspect  ratios  and  the  distortion  factors  of  the  elements. 

t  NASTRAN  MSGMESH  is  developed  by  MacNeal-Schwendler  Corporation 
□  GIFTS  is  developed  by  Sperry  Univac  Computer  System 
•  SUPERSAP  is  developed  by  Algor  Corporation 
f  I-DEAS  is  developed  by  Structural  Dynamic  Research  Corporation. 
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Although  these  pre-processors  help  in  generating  acceptable  meshes, 
it  is  still  difficult  to  obtain  a  mesh  that  is  best  suited  for  the  problem  at 
hand.  The  difficulty  lies  in  the  definition  of  the  "  best  "  mesh  .  Is  there 
a  best  mesh  for  a  particular  domain  ?  If  so,  is  there  a  different  one  for 
different  set  of  boundary  conditions  or  a  different  set  of  loading  ?  Is 
there  a  different  optimal  mesh  for  different  differential  equations  in  the 
same  domain  ?  Answers  to  these  questions  are  discussed  in  the  following 
sections. 


1.3  OPTIMAL  MESH 

Recall  from  section  1.1  that  the  functional  is  a  function  of  the 
unknown  or  dependent  variables.  Note  that  it  is  also  a  function  of  the 
coordinates  of  the  nodes.  Therefore  it  can  be  expressed  as: 

(1.2)  t r  =  7r(ui,dk) 

where,  Uj  is  the  vector  of  unknown,  dk  is  the  position  vector  of  kth 
node. 


In  order  to  obtain  a  true  minimum  on  (1.2),  in  addition  to  the 
equilibrium  equations  (1.1),  it  is  necessary  to  consider  the  following 
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equations. 


(13) 

frr  _  1  dKjj 

ddk  2  Ui  adk  Uj 


where,  rk  is  the  residual  vector. 


Solution  of  (1.3)  along  with  the  geometrical  constraints  will  yield  the 
optimal  locations  of  the  nodes,  which  when  used  in  (1.1)  should  result  in 
a  uj  that  is  closest  to  u^cf 


The  method  seems  to  be  very  simple,  theoretically.  However,  the 
non-linear  algebraic  equations  (1.3)  are  difficult  to  solve  explicitly.  Even 
for  a  simple  geometry  in  one  dimension  the  algebra  is  very  complicated. 
Numerical  solutions  are  also  difficult  [8].  Some  of  the  solution  methods 
for  non-linear  equations  like  gradient  methods  and  complex  methods  have 
been  tried  but  with  little  success.  Among  investigators  examining  this 
problem,  Prager[9]  has  made  a  note  worthy  contribution.  He  examined  a 
bar  with  a  linearly  varying  cross  section  under  tension.  He  showed  that 
the  grid  producing  the  desired  least  potential  energy  is  the  one  where  the 
cross  section  areas  at  the  nodes  form  a  geometric  series.  This  problem  is 
studied  in  greater  detail  in  the  next  chapter. 
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1.4  MESH  REFINEMENT 


As  described  in  section  1.2  the  user  needs  to  select  the  number  of 
nodes  and  elements  in  the  model.  The  selection  may  be  the  one  that 
leads  to  the  best  description  of  the  domain  geometrically.  For  example,  a 
curved  surface  could  be  modelled  by  a  series  of  interconnected  flat 
rectangular  facets.  The  larger  the  number  of  facets,  the  better  is  the 
model.  The  selection  may  also  be  based  upon  intuition,  past  experience 
and  engineering  judgement.  The  mesh  obtained  may  be  adequate  in  some 
cases.  In  other  cases,  especially  when  singularities  are  present,  the  mesh 
may  not  be  adequate  to  obtain  the  results  to  the  accuracy  desired.  In 
such  cases,  the  meshes  need  to  be  refined. 


1.4.1  REFINEMENT  PROCESS 


There  are  three  ways  of  refining  a  finite  element  mesh: 

a)  The  H-method:  This  method  increases  the  number  of  elements 
and  hence  decreases  the  element  size  while  keeping  the  polynomial  order 
of  the  shape  function  constant. 

b)  The  P-method:  This  method  increases  the  polynomial  order  of 
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the  interpolation  function  while  keeping  the  number  of  elements  in  the 
model  constant. 

c)  The  R-method:  This  method  redistributes  the  nodes  while  keeping 
the  element  number  and  the  polynomial  order  of  the  interpolation 
function  constant. 


1.4.1. 1  H  -  Method 

This  method  is  primarily  based  on  the  choice  of  characteristic  length 
of  the  elements.  "Characteristic  length  "  is  referred  to  in  a  generalized 
sense  and  is  required  to  define  the  element  topologically.  A  linear 
element  requires  one  characteristic  length,  whereas  an  element  of 
rectangular  shape  requires  two  characteristic  lengths  and  a  triangular 
element  requires  three  characteristic  lengths  for  its  definition.  In  the 
triangular  element  the  three  length  informations  may  be  any  combination 
of  lengths  and  angles. 

Instead  of  expressing  the  functional  in  terms  of  the  position  vectors 
of  the  nodes,  as  in  (1.2),  it  can  be  expressed  as  a  function  of  the 
element  characteristic  lengths  as  : 
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(1.4) 


7T  =  TrCUj.hik) 


where,  is  the  element  characteristic  length,  1  is  the  index  on  the 

characteristic  length  for  element  k 

Also,  note  that  there  will  be  geometrical  constraints  on  h^.  For 
example,  the  sum  of  the  element  lengths  in  a  particular  direction  should 
be  equal  to  the  overall  dimension  of  the  model  in  that  direction.  Again 
as  described  in  section  1.3,  the  function  can  be  minimized  with  respect  to 
the  characteristic  lengths. 


(1.5) 


(h r 
dhlk 


1 

2Ui 


3Kij 

dhik 


di 

dhlk 


uf  = 


rk  =  0 


Solving  (1.5)  along  with  the  constraints  yield  the  characteristic 
lengths  and  hence  defines  the  best  mesh.  Equation  (1.5)  is  equivalent  to 
(1.3)  cast  in  the  frame  work  of  characteristic  lengths.  Therefore  the 
solution  as  indicated  in  section  1.3  is  difficult.  A  practical  procedure 
using  this  method  consists  of  selecting  a  coarse  initial  mesh,  solving  the 
equilibrium  equations  and  computing  the  residue  rk  on  each  element.  The 
set  of  elements  with  large  values  of  residues  is  the  region  that  needs  to 
be  refined.  The  identified  region  can  be  refined  by  sub-dividing  the 
elements,  thus  creating  new  regions,  or  by  deleting  all  the  elements  in  the 
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region  and  replacing  them  by  finer  elements.  However,  the  new  elements 
need  to  be  of  the  same  type  as  those  in  the  initial  mesh.  The  equations 
of  the  new  model  are  solved  and  the  residues  are  computed.  If  the  values 
of  the  residues  are  still  large,  the  refinement  procedure  can  be  repeated. 
Indeed  ,  it  could  be  used  iteratively  until  the  solution  meets  the 
prescribed  accuracy. 

The  monotonic  convergence  of  the  refinement  procedure  has  been 
studied  by  Melosh  [10]  and  Key  [11].  A  convergence  theorem  has  been 
introduced  by  Carroll  and  Baker  [12],  which  states  : 

Theorem:  A  necessary  consequence  of  the  following  refinement  sequence 

(1.6)  7Tn  >  7Tn>1  >  >  •  ’  •  > 

where,  m  represents  successive  refinements  of  the  initial  finite  element 
mesh  n,  is  the  existence  of  an  optimum  sub-division  such  that 

(1.7)  (hi)  <  (hik) 

where  hi  correspond  to  the  optimum  mesh. 

The  usefulness  of  this  theorem  can  be  explored  in  the  discussion  of 
the  r  -  method.  The  difficulty  in  using  this  method  is  in  the  estimation  of 
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the  derivatives  involved  in  the  computations  of  the  residues. 


1.4. 1.2  P  -  Method 

This  method  is  primarily  based  on  the  choice  of  the  order  of  the 
interpolation  function,  which  in  practice,  translates  to  the  choice  of 
element  type.  For  example,  in  a  two  dimensional  domain,  the  basic 
triangular  element  with  three  nodes  at  the  three  vertices  uses  a  linear 
shape  function  (p=l).  In  order  to  choose  quadratic  shape  functions  (p=2), 
the  triangular  element  with  six  nodes,  three  at  the  vertices  and  three  at 
mid-side  locations,  has  to  be  selected.  Similarly,  for  cubic  functions,  an 
element  with  nine  nodes  is  selected. 

Higher  order  elements  generally  provide  better  description  of  the 
domain  geometrically.  They  are  particularly  useful  in  regions  where  the  use 
of  lower  order  elements  would  result  in  a  mesh  with  poor  aspect  ratios  in 
those  elements.  From  the  point  of  view  of  solution  accuracy,  higher  order 
elements  are  usually  more  accurate  than  the  lower  order  elements.  But 
this  does  not  mean  that  increasing  the  polynomial  order  indiscriminately 
will  always  provide  point-wise  convergence  to  the  exact  solution.  The 
argument  is  based  on  the  theory  of  interpolation.  Prenter  [13]  states  that 
this  notion  on  convergence  was  first  dispelled  by  Meray  and  later  by 
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Runge.  He  illustrates  this  with  the  function  f(x)  =  l/(  1+5x2)  being 
interpolated  by  Lagrange  polynomial  of  order  5  and  15  with  evenly  spaced 
nodes  in  the  interval  [-1,1]  which  display  divergence  at  -  1  and  1. 
Although  the  example  is  for  a  continuous  interpolation  function  rather 
than  a  piecewise  function,  as  in  a  finite  element  model,  it  shows  that 
there  is  good  reason  to  exercise  caution  in  increasing  the  polynomial 
order. 


1.4.1.3  R  -  Method 

This  is  a  far  less  explored  method.  It  neither  increases  the 
polynomial  order  nor  decreases  the  element  character  length.  The  mesh  is 
refined  simply  by  re-distributing  the  nodes  in  the  domain  such  that  the 
potential  energy  is  reduced. 

Recall  the  theorem  introduced  by  Carroll  and  Baker,  stated  in 
section  1.4.1. 1.  The  theorem  indicates  that  there  exists  an  optimal 
distribution  of  the  nodes  in  a  domain.  Any  other  distribution  will  yield  a 
potential  energy  higher  than  the  lowest  possible  for  the  given  number  of 
degree  of  freedom.  The  theorem  also  indicates  that  :  given  a  distribution 
of  nodes,  a  new  distribution  will  be  a  refinement  over  the  old  distribution 
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only  if  it  results  in  a  lower  potential  energy  than  the  old  distribution. 
This  fact  could  be  used  in  an  iterative  refinement  procedures. 


1.4.2  ADAPTIVE  MESH  REFINEMENT 


The  refinement  that  follows  the  requirements  of  the  differential 
equation  or  the  boundary  conditions  closely  is  called  an  adaptive 
refinement.  This  method  is  used  to  tailor  the  mesh,  including  finer 
elements.  It  can  also  provide  elements  of  higher  polynomial  order  where 
necessary  as  opposed  to  the  method  of  h  or  p  *  refinement  all  over  the 
domain.  The  practical  method  mentioned  in  section  1.4.1. 1  is  adaptive.  The 
obvious  advantage  is  that  it  achieves  the  desired  accuracy  level  while 
keeping  the  number  of  degrees  of  freedom  low. 


1.4.3  A  -  PRIORI  AND  A  -  POSTERIORI  METHODS 


The  classification  of  methods  into  a-priori  and  a-posteriori  refers  to 
refinement  before  and  after  the  solution  of  the  equilibrium  equations.  In 
a  finite  element  program  the  solution  process  is  one  that  needs  much  of 
computer  time.  If  discretization  errors  can  be  estimated  a-priori,  then  the 
mesh  can  be  suitably  altered  to  obtain  the  best  accuracy  possible  by 
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solving  the  equations  only  once.  Unfortunately  there  are  no  practical  a* 
priori  methods  available.  The  author  has  not  found  any  in  the  literature 
survey.  This  study  is  an  attempt  to  provide  one.  There  are  several  a- 
posteriori  methods  available  for  refinement. 


1.4.4  USE  OF  HIERARCHICAL  ELEMENTS  IN  REFINEMENT 


An  hierarchical  displacement  element  is  one  whose  stiffness  matrix 
contains  the  stiffness  matrices  of  lower  order  elements  explicitly  as 
submatrices[14]. 

Consider  a  two-node  axial  element.  Its  stiffness  matrix  is  given  by: 

A  E  rl  -1 

1  l-l  1 


An  hierarchical  displacement  element  of  one  higher  order  contains 
an  additional  node  in  the  middle  of  the  element,  as  in  the  conventional 
quadratic  element.  However  the  shape  function  chosen  for  the  midside 
node  is  different  from  the  one  chosen  in  the  conventional  element.  This 
results  in  the  stiffness  matrix: 


A  E 
1 


1 

-1 

0 


-1 

1 

0 


0 

0 
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Note  that  the  stiffness  matrix  of  the  basic  element  is  contained  in 
the  new  matrix  as  a  submatrix.  The  stiffness  matrices  of  higher  order 
elements  are  built  by  a  similar  process  if  a  higher  order  element  is  coded 
into  the  finite  element  program,  it  includes  stiffness  matrices  of  all  lower 
order  elements.  In  the  process  of  refinement  if  a  higher  order  element  is 
chosen,  the  previously  computed  stiffness  coefficients  would  still  be  valid. 
Hence,  only  a  few  additional  coefficients  have  to  be  evaluated.  The 
method  is  easier  than  the  conventional  p  -  method  of  increasing  the 
polynomial  order  where  the  computation  of  the  entire  higher  order 
element  stiffness  matrix  is  required. 

Refinement  using  hierarchical  elements  is  a-posteriori  and  appears  to 
be  attractive.  However  more  research  work  needs  to  be  done  in  this  area. 

1.5  RESEARCH  OBJECTIVES 

It  is  clear  that  the  accuracy  of  the  finite  element  results  is  mesh 
dependent.  A  proper  mesh  selection  procedure  is  therefore  necessary.  A 
posteriori  methods  are  adaptive  in  nature  but  are  expensive  in  terms  of 
computer  processing  time.  On  the  other  hand,  a  priori  methods  are  not 
adaptive.  They  use  geometrical  criterions,  element  aspect  ratios,  for 


15 


example,  for  improvements.  Some  of  them  help  estimating  the  overall  error. 
They,  however,  do  not  indicate  the  regions  which  need  refinement.  The 
prime  objective  of  this  report  is,  therefore,  to  develop  a  criterion  which  helps 
in  identifying  the  region  for  refinement  or  rezoning  process  even  before  the 
equilibrium  equations  are  solved.  The  procedure  based  on  the  criterion 
should  be  able  to  guide  the  user  in  improving  the  grids. 

Finally,  the  report  itself  is  based  upon  the  first  author’s  doctoral 
dissertation  [15]  at  the  University  of  Cincinnati. 
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2.  ANALYSIS 


2.1  ANALYTICAL  APPROACH 


The  objective  is  to  develop  a  practical  and  efficient  procedure  of 
grid  enhancement  and  optimization.  The  thesis  is  that  for  many  problems 
the  minimization  of  the  trace  of  the  stiffness  matrix  with  respect  to  the 
nodal  coordinates,  leads  to  a  minimization  of  the  potential  energy,  and  as 
a  consequence,  provides  the  optimal  grid  configuration.  To  see  this, 
consider  the  governing  matrix  equation  of  finite  element  analysis: 

(2.1)  [K]{u}  =  {[} 

where,  [K]  is  the  stiffness  matrix,  {u}  is  the  array  of  dependent 
variables,  and  {f}  is  the  force  array. 

Matrix  [K]  can  be  viewed  as  an  operator  which  maps  {u}  into  {f}.  In  this 
context,  since  [K]  is  symmetric,  an  orthogonal  transformation  {T},  which 
diagonalizes  [K],  can  be  found.  That  is, 

(2.2)  [K]  .  [T]t[K][T] 
where  [K]  is  a  diagonal  matrix. 
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Let  [T]{u}  and  [T){f}  be  {u}  and  {?}.  Then  the  potential  energy  n 
may  be  expressed  as: 

(23)  *  =  j  {u}t[KHu}  -  {f}T{u} 

In  terms  of  the  array  components,  r  becomes: 

(2.4)  *  -  t  4  kafif  '  4«i 

*=K  ^ 

A  A 

where  the  (i=l,2,...,n)  are  the  diagonal  entries  of  [K] 

n  ~  ^ 

Observe  in  Equation  (2.4)  that  the  last  term:  J^f^  does  not  explicitly 

i=l 

involve  the  nodal  coodinates.  Therefore,  it  does  not  effect  the 
minimization  of  x  with  respect  to  the  nodal  coordinates.  Also,  since  the 
u2  are  positive  and  are  independent  variables  in  the  minimization  of  n, 
the  minimization  of  n  with  respect  to  the  nodal  coordinates  occurs  when 
the  sum  of  the  ky  (the  trace  of  [k])  is  a  minimum.  Since  the  trace  of  a 
matrix  is  invariant  an  orthogonal  transformation,  minimizing  the 

trace  of  [K]  is  equivalent  to  minimizing  the  trace  of  [K]. 
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2.2  ENERGY  APPROACH 


Consider  a  one  -  degree  of  freedom  system.  The  external  work  done 
(=fu)  varies  linearly  with  respect  to  u.  Also  the  strain  energy  (  =  1/2  Ku2) 
varies  quadratically  with  respect  to  u.  Potential  energy  is  the  difference  of 
strain  energy  and  work  done.  See  Figure  2.1. 

From  the  instant,  the  structure  is  loaded  the  operating  point  moves 
from  the  origin  to  the  point  where  the  potential  energy  reaches  its 
minima  (equilibrium). 

Now  consider  the  structure  with  a  reduced  stiffness  (K’<  K).  The 
new  strain  energy  and  the  potential  energy  curves  are  also  shown  in  the 
figure.  Note  that  the  strain  energy  curve  has  become  slightly  flatter. 
Therefore  the  potential  energy  curve  has  reached  a  new  low,  which  is 
lower  than  the  previous.  The  displacement  has  improved  from  ua  to  Ub- 
Therefore  it  is  quite  clear  that  in  a  single  degree  of  freedom  case,  a  less 
stiff  structure  has  a  lower  potential  energy  than  the  stiff  one. 

Next,  consider  an  n  -  degree  of  freedom  system.  Using  an 
orthogonal  transformation  matrix,  [K]  can  be  diagonalized.  This  would  de¬ 
couple  the  degrees  of  freedom.  Therefore  each  degree  of  freedom  can  be 
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compared  with  the  single  degree  of  freedom  system  as  described  above.  If 

A  A 

[K]  is  the  transformed  stiffness  matrix,  then  finding  minimum  ku 
(i=l,2,...n)  would  yield  the  best  grid.  Since  for  the  mesh  configuration  the 

A  A  A 

minimum  k,j  have  been  found,  the  trace  of  [K]  which  is  the  sum  of  k^ 
will  also  be  a  minimum.  But  the  trace  is  an  invariant  under  orthogonal 
trasformation.  Therefore  minimization  of  trace  of  [K]  would  lead  to 
minimization  of  potential  energy. 

It  should  be  noted  that  the  diagonal  dominance  of  [K]  is  not 
adversely  affected  by  the  minimization  of  trace.  The  improved  stiffness 
matrix  is  the  result  of  redistribution  of  the  nodes  and  of  not  any  arbitrary 
mathematical  operation. 
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3.  APPLICATIONS 


Applications  of  concepts  developed  are  illustrated  in  the  examples  in 
this  chapter. 


3.1  TAPERED  BAR 


3.1.1  DESCRIPTION 


Consider  the  axially  loaded  tapered  bar  shown  in  Figure  3.1.  (This  is 
the  same  problem  examined  by  Prager  [9]  and  Masur  [16].  The  objective  is 
to  determine  a  finite  element  mesh  which  best  predicts  the  axial 
displacement.  Let  the  bar  have  length  L  and  let  it  be  divided  into  n  elements 
with  n  + 1  nodes  (numbered  0  to  n)  as  shown.  Let  the  areas  at  the  ends  of 
the  bar  be  A„  and  A[.  Let  £  be  the  non  dimensional  length  parameter 
defined  as: 

(3.1)  x 

e  =  - 

L 

Then  the  area  at  any  particular  £  along  the  bar  is 


22 


(3.2) 


A  =  Ao(l-ce) 


where  c  is 

(3.3) 


c 


Aq-Ai 

Aq 


Hence,  the  area  at  the  ith  node  is: 

(3.4)  Ai  =  A0(Ki) 


where,  &  is  X;/L. 

Let  the  individual  elements  have  a  uniform  cross  section.  For 
example,  let  the  ith  element  have  cross  section  area  A;  and  length  1;  as  in 
Figure  3.2.  (Note  that  the  elements  do  not  necessarily  have  the  same 
length.)  Then  A;  and  lj  are: 


—  Aj.i+Ai 

AS  =  — 

(3.5) 


1;  =  Xj  Xj.x  -  L(£j  "  £i-i) 


The  element  stiffness  matrix  for  the  ith  element  is  [5]  : 


(3.6) 


M 


AiE  r  ! 

li  l'1 


*1 

1 


] 


where  E  is  the  elastic  modulus.  Then  the  trace  r  of  the  global  stiffness 
matrix  is: 


(3.7) 


T 


2E  £ 


»=i 


Ai 


'Ah  ♦Af 

.6  -  e,i. 


The  improved  node  location  occurs  when  the  trace  is  minimized 
with  respect  to  the  nodal  coordinates  £i  (i=l,2,...n-l).  Hence,  by  setting  the 
derivative  of  r  with  respect  to  &  equal  to  zero  we  obtain: 


(3.8) 


parameter  r^  defined  as  lj/lj.  Then  the  ratio  may  be  written  as: 


(3.10) 


i-t-i  U 

h  =  Ji 

ii 


ri+l,l 

til 


Then  L/lj  is 


(3.11) 


t.Ai.f  !l.!i 

1.  "  ft  li  U  <n  r„ 


where  Sn  is  defined  as  ^r^. 

H 


Using  this  notation  Equation  (3.9)  may  be  rewritten  as 


(3.12) 


Aj+i  =  A; 


rM,l 


Til 


+  Aj.j 


ri+l,l 


til 


cApr^i 

S„ 


ri+l,l 


til 


+  1 


Also,  Ci  may  be  written  as: 


(3.13) 


6  = 


Ii  +  h 


25 


1  +  r21  +  r3i  + 


1  ‘  _  Sj_ 

S  O  S  “  C 

sn  jci  an 


Then,  from  equation  (3.4),  Aj  may  be  written  as 


(3.14) 


Aj  *  A0  1-c  — 


Specifically,  Ax  and  A2  are 


A,*Ao  1  - 


(3.15) 


A2  =  Aq  1  -  c 


(1  +  r21) 


To  obtain  the  element  area  ratios  let  i=l  in  equation  (3.12).  A2  is 


(3.16) 


A2  =  A^r21-  1  j  +  Aor2i  +  cAor21  - 
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Then,  by  using  equations  (3-15)  to  eliminate  Ax  and  A2,  we  obtain 


(3.17)  Sn(l  -  r21)  =  c 


From  the  first  of  equations  (3.15)  we  have 


(3.18) 


Aj  =  Aor2J 


or 


Similarly,  the  second  of  equations  (3.15)  leads  to 


A2  =  A0r2i 


(3.19) 


or 


A2 

Ai 


=  r21 


Next,  let  i=2  in  equation  (3.12).  By  using  the  same  procedure  we 
obtain 


(3.20) 


1  '  r32  =  (1  *  r32  +  r2I 


) 
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But,  in  view  of  equation  (3.17)  this  becomes 

1  *  *32  =  0  '  r2l)  (1  '  r32  +  r2l) 

(3.21)  or 

r2l(r32  '  r2l)  s  0 


Thus, 

(3.22)  r32  a  r2j 

From  equation  (3.14)  we  have 


Therefore,  we  obtain 


(3.24) 


A3  Aj 
A2  Aj 


Proceeding  similarly  for  i=3,4,...  we  get 
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(3.25) 


An  An.j 

A.n-1  An.2 


Ai 

— —  =  r21  =  r32  =  •  •  •  =  rn>n.j 


Thus,  we  have  the  relations 


2 

r31  =  r21r32  =  r21 


(3.26) 


3 

r41  =  r43r32r21  =  r2l 


fil 


ri-l 

r21 


Hence,  S;  is  the  geometric  series 


(3.27) 


Sj  =  1  +  r21  +  r21  + 


ri-i  _ 
r21  - 


1  •  r21 
(1  ‘  r2l) 


Then,  from  equations  (3.14)  and  (3.17),  we  have 


(3.28) 


A;  =  Aq^I  -  (1  -  r2l)  |  =  A0r2i 

and  then 


A  n  =  A0r2i  =  A; 
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Then,  from  equation  (33)  we  see  that  t2l  is 


(3.29)  r„  -  (1  -  c)>/« 


Finally,  by  substituting  into  equation  (3.13)  we  have 


(330) 


=  Y  (1  +  r*i  +  131  +  '  '  '  +  ri») 


=  (1  *  rai) 


f  1  +  r2i  +  th  + 


*  ri-l'l 
+  r2i 


1  •  fai  _  1  -  (1  -  c)fr 
c  c 


This  is  the  result  obtained  by  Prager  [9]  in  his  analysis  of  the  same 
problem. 


3.1.2  DISCUSSION 


First,  observe  that  in  Equation  (330)  for  a  uniform  thickness  beam 
c  is  zero  and  thus  is  undetermined.  This  means  that  for  a  uniform 
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thickness  beam  the  nodal  position  are  arbitrary.  That  is,  all  mesh  are 
equally  optimal  for  a  uniform  thickness  beam. 

Next,  consider  again  the  element  stiffness  matrix  of  equation  (3.6). 
From  equations  (3.4)  and  (330)  the  scalar  multiplier  is 


(3.31) 


EAj  E(Ai.!  +  Aj)  AqC  f  1  +  r2i 

T~ =  2L(e;  -  eM)  =  IT  i  -  r2T 

/ 


Since  this  is  a  constant  (independent  of  i)  the  element  stiffness  matrix 
is  the  same  for  each  element.  This  means  that  each  element  has  the  same 
strain  energy.  Masur  [15]  has  suggested  that  this  result  is  due  to  the  simple 
geometry  of  the  problem. 


Even  with  this  simple  geometry  however,  like  the  methods  discussed 
in  section  13,  the  analysis  needed  to  determine  the  optimal  nodal 
positions  has  been  extremely  detailed.  With  more  complex  geometries  the 
analysis  will  become  intractable.  However,  it  is  not  necessary  to  obtain 
recursive  relationships  by  analytical  methods  as  employed  in  this  example. 
The  criteria  of  minimizing  the  trace  of  the  stiffness  matrix  is  a 
comparatively  simple  procedure  -  readily  amenable  to  the  development  of 
computer  algorithms  for  optimal  nodal  locations. 
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3.1.3  NUMERICAL  EXAMPLE 

To  illustrate  the  value  of  optimizing  the  mesh  consider  an  axially 
loaded  bar  which  tapers  to  1/3  the  base  area  as  in  Figures  3.1  to  3.3. 
Specifically,  let  P,  A*,  C,  E,  and  L  have  the  values: 

(3.31)  P  =  20N,  A.  =  0.0015m:,  C  =  2/3 
E  =  2.07  x  10“N/mJ,  L  =  4m 


The  objective  is  to  find  the  axial  displacement.  From  elementary 
mechanics  the  axial  displacement  u  at  any  location  x  is: 


(3.33) 


u 


PL 

AoEc 


In 


1- 


£x 

L 


To  compare  the  displacement  results  of  finite  element  models  with 
Equation  (3.33),  four  models  of  the  bar,  each  having  four  elements,  were 
examined.  One  of  the  models  had  a  uniform  nodal  distribution.  Another  had 
the  "optimal"  mesh  as  developed  in  Equation  (3.33).  The  remaining  two 
models  had  arbitrarily  selected  nodal  distributions.  The  nodal 
displacements  were  evaluated  using  the  four  models  and  compared  with 
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the  displacements  calculated  by  (333).  Table  3.1  shows  the  results.  Table 
3.2  presents  an  error  analysis  and  also  an  L?  -  norm  of  the  error.  As 
expected  the  optimum  mesh  produces  the  least  Lj  error. 


TABLE  3.1  -  Comparison  Of  Axial  Displacements  For  The  Tapered  Bar  Of 
Figure  3.1  Calculated  Using  Various  Models 


Axial 

location 

x,  m 

Exact 

displacement 
eq.  (3.30), 
10'9m 

Displacements  using  various  models,  10'9m 

Uniform 

Mesh 

"Optimum" 

Mesh 

(Prager) 

Mesh  3 

Mesh  4 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.5 

33.6276 

33.60639 

33.60639 

1.0 

70.46244 

70.2679 

70.41338 

1.4409860 

106.1461 

105.4820 

2.0 

156.7020 

156.1509 

15.56506 

2.5 

208.3078 

206.8158 

207.1804 

2.5358986 

212.2923 

210.9648 

3.0 

267.8830 

266.5719 

3.3678522 

318.4384 

316.4529 

4.0 

424.5845 

421.1612 

421.940 

417.6195 

417.9814 

TABLE  3.2  -  Error  Analysis  (Tapered  Bar) 
Error  for  various  meshes,  10'9m 


Axial 

1  Uniform 

"Optimum" 

Mesh  3 

Mesh  4 

location 

mesh 

mesh 

x,  m 

11K2I2H 

0.0 

0.5 

0.0 

0.0 

0.0 

0.02121 

1.0 

1.4409860 

0.1945411 

0.664118 

0.04906 

2.0 

0.5506105 

1.0514 

2.5 

2.5358986 

1.32769 

1.492 

1.1274 

3.0 

3.3678522 

1311108 

1.985439 

4.0 

3.423228 

2.64433 

6.965 

6.6031 

12  -  norm 

3.7119418 

3.6246009 

7.1232118 

6.780697 

3.1.4  OBSERVATIONS 


[1]  The  analysis  and  the  numerical  results  demonstrate  the  potential 
usefulness  of  the  trace  minimization  mesh  improvement  method. 

[2]  Minimization  of  the  trace  of  the  stiffness  matrix  is  a  relatively 
simple  mesh  optimization  procedure.  It  is  readily  adaptable  to 
algorithm  development. 
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Figure  3.3  *  Tapered  Bar  with  Axial  Load 
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3.2  HEAT  TRANSFER  IN  AN  INFINITE  CYLINDER 


3.2.1  CONFIGURATION  AND  PROBLEM  DEFINITION 


Consider  an  annular  cylinder  with  Infinite  length  having  inner  and 
outer  radii:  r0  and  rn.  Let  the  thermal  conductivity  be  /c.  Let  the 
temperatures  at  the  inner  and  outer  radii  be:  T0  and  Tn.  Then  the 
governing  equation  for  the  temperature  distribution  along  a  radial  lines  is: 


(3.34) 


d  dt  n 

—  r/c—  =  0 

dr  dr 


The  boundary  conditions  are: 

(3.35)  T  =  T0  at  r  =  r0 

T  =  Tn  at  r  =  r„ 


The  solution  of  equation  334  subject  to  equation  3.35  is: 


(336) 


T  =  Tn  +  (T0  * 
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Next,  suppose  that  the  temperature  gradient  at  the  inner  surface  is 


specified  as:  q0.  The  boundary  conditions  are  then: 


(3.37) 


at  r  =  r0 


T  =  Tn  at  r  =  rn 


In  this  case  the  solution  of  equation  334  is: 


(338) 


T  =  Tn 


r0qoln 


3.2.2  FINITE  ELEMENT  FORMULATION  AND  MESH  OPTIMIZATION 


Figure  3.4  shows  the  finite  element  model.  It  consists  of  a  series  of 
annular  elements.  For  elements  (e)  let  the  inner  and  outer  radii  be  re 
and  r^.  The  entries  of  the  stiffness  matrix  are: 


(3.39) 


kjj  =  2/kk^  J* 

r. 


'  dN*' 

f  dN®' 

dr 

> 

dr 
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where  /c,  are  the  element  conductivity  constants  and  where  the 


element  shape  functions  Nf  and  N*  are: 


(3.40) 


and 


(r  •  r‘) 
(re+l  *  r«) 


By  carrying  out  the  indicated  operations  the  element  stiffness  matrix 
becomes: 


(3.41)  [  k}  ]  -  S.  [,\  V  ] 


where  Se  is  defined  as: 


(3.42) 


Se 


7T/Ce 


U  *  Te-1 
ee  *  re-l 


Hence  the  trace  r,  of  the  global  stiffness  matrix  is: 


(3.43) 


where  n  is  the  the  number  of  elements. 
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The  trace  may  be  minimized  with  respect  to  the  nodal  coordinates 
by  setting  the  partial  derivative  of  r  with  respect  to  re,  equal  to  zero. 
This  leads  to  the  relatively  simple  relation: 


(3.44) 


-  -L- 

r.  "  re-l 


By  repeated  use  of  this  relation  the  nodal  positions  are  given  by: 


(3.45) 


re  -  *o 


e/n 


3.2.3  NUMERICAL  EXAMPLES 

To  illustrate  the  effectiveness  of  the  method  considers  the  annular 
cylinder  with  the  following  temperatures  specified  on  the  boundaries: 

(3.46)  r0  =  20  mm  T0  =  100°  C 

rn  =  50  mm  Tn  =  0°  C 
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Let  the  conductivity  be  constant  throughout  the  cylinder  with  value:  1.0. 

Consider  two  finite  elements  models,  each  with  four  elements:  Let 
the  first  have  a  uniformly  spaced  mesh.  Let  the  second  have  a  mesh  with 
nodal  spacing  governed  by  equation  3.45.  The  objective  is  to  determine 
the  temperature  distribution  across  the  thickness. 

The  solution  of  the  finite  element  governing  equations  lead  to  the 
results  listed  in  Table  3.3.  (The  temperatures  at  the  intermediate  points,  if 
they  are  not  obtained  directly,  are  obtained  using  linear  interpolation 
between  the  nodal  values.)  The  error  is  defined  as  the  difference  between 
the  theoretical  results  and  the  finite  element  results.  The  mesh  governed 
by  equation  334  (called  the  "improved"  mesh)  is  found  to  have  zero 
errors  at  the  nodes.  Hence,  the  norm  of  the  errors*  is  much  smaller 
than  that  of  the  uniform  mesh. 
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TABLE  3.3  -  Comparison  Between  Finite  Element  Uniform  Mesh  and  Improved 
Mesh  Temperature  Results  with  Theoretical  Values  for 
Temperature  Specified  Boundary  Conditions 


Radius 

(mm) 

Temperature 

C 

Error 

Uniform 

Mesh 

Improved 

Mesh 

Theoretical 

Values 

Uniform 

Mesh 

Improved 

Mesh 

20.0 

100.0 

100.0 

100.0 

0.0 

0.0 

25.148669 

76.21655 

75.0 

75.0 

-1.216559 

0.0 

27.5 

65.354967 

65.920251 

65.24534 

-0.109625 

-0.674909 

31.622777 

50.881148 

50.0 

50.0 

-0.881148 

0.0 

35.0 

39.024743 

39.628662 

38.92596 

-0.098784 

-0.702703 

39.763536 

25.538185 

25.0 

25.0 

-0.538185 

0.0 

42.5 

17.790692 

18.316874 

17.73661 

-0.05408 

-0.580262 

50.0 

0.0 

0.0 

0.0 

0.0 

0.0 

Lj  norm  of  error  : 

1.6033656 

1.1340184 

Next,  consider  the  same  cylinder  but  let  the  temperature  gradient  be 
specified  on  the  inner  boundary.  Specifically,  let  the  boundary  conditions 
be: 

(3.46)  4^  =  -5.4567833’  C/mm  at  r0  =  50  mm 

dr 

T  =  T„  =  O’C  at  rn  =  50  mm 

(The  temperature  gradient  of  -5.4567833  ’  C/mm  on  the  inner  boundary 
leads  to  the  same  theoretical  temperature  distribution  as  in  the  first 
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example.)  Table  3.4  shows  the  comparisons  between  the  finite  element 
solutions  and  the  theoretical  values  of  the  temperature.  Once  again  the 
norm  of  the  error  shows  that  the  improved  mesh  provides  better 
results  than  the  uniform  mesh. 


TABLE  3.4  -  Comparison  Between  Finite  Element  Uniform  Mesh  and  Improved 
Mesh  Temperature  Results  with  Theoretical  Values  for 
Temperature  /  Temperature  Gradient  Specified  Boundary 
Conditions 


Temperature  '  C  Error 


Radius 

(mm) 

Uniform 

Mesh 

i - 

Improved 

Mesh 

Theoretical 

Values 

' - 

Uniform 

Mesh 

Improved 

Mesh 

20.0 

99.47716 

99.570031 

100.0 

0.52284 

0.42997 

25-148669 

75.818072 

74.677527 

75.0 

-0.818072 

0.322473 

27.5 

65.01327 

65.636818 

65.24534 

0.23207 

-0.391478 

31.622777 

50.615125 

49.78502 

50.0 

-0.615125 

0.214981 

35.0 

38.82071 

39.458273 

38.92596 

0.10525 

-0.532313 

39.763536 

25.404668 

24.892511 

25.0 

-0.404668 

0.10749 

42.5 

17.69768 

18.238117 

17.73661 

0.03893 

-0.501507 

50.0 

0.0 

0.0 

0.0 

0.0 

0.0 

Lj  norm  of  error : 

1.245467 

1.0172293 

Since  the  temperature  gradient  is  specified  at  the  inner  boundary,  it 
is  also  of  interest  to  know  how  the  values  of  the  temperature  gradients 
obtained  using  the  two  finite  element  meshes  compare  with  each  other 


and  with  the  theoretical  values.  Table  3.5  provides  such  a  comparison, 
(the  temperature  gradients  at  the  intermediate  points,  if  they  are  not 
obtained  directly,  are  computed  by  using  forward  differences.)  The 
improved  mesh  has  a  small  error  at  the  inner  radius  as  well  as  a  smaller 
L2  norm  of  errors  overall. 

3.2.4  DISCUSSION 


The  numerical  example  show  that  the  "improved"  mesh  of  equation 
3.45  produces  results  which  are  closer  to  the  theoretical  values  than  those 
obtained  using  the  uniform  mesh.  Therefore,  the  mesh  of  equation  3.45  is 
an  improvement  over  the  uniform  mesh  for  both  the  temperature  fixed 
boundary  conditions  and  for  the  mixed  boundary  conditions. 

The  values  of  stiffness  matrix  traces  of  the  uniform  and  improved 
meshes  are  74.666666  and  70.151974  respectively.  This  indicates  that  for 
these  examples  the  trace  is  not  especially  sensitive  to  the  nodal  locations. 
Therefore,  the  difference  in  the  Lj  norms  of  error  are  not  great.  More 
dramatic  difference  in  the  results  will  occur  in  problems  where  the  trace  is 
more  sensitive  to  changes  in  the  nodal  positions. 
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TABLE  3.5  -  Comparison  Between  Finite  Element  Uniform  Mesh  and  Improved 
Mesh  Temperature  Gradient  Results  with  Theoretical  Values  for 
Temperature  /  Temperature  Gradient  Specified  Boundary 
Conditions 


Radius 

(mm) 

Temperature  *  < 

Error 

Uniform 

Mesh 

Improved 

Mesh 

Theoretical 

Values 

Uniform 

Mesh 

Improved 

Mesh 

20.0 

-45951853 

■ 

-4.8347453 

•5.4567833 

0.861598 

0.622038 

25.148669 

•4.5951853 

-3.8449325 

-4.3396146 

-0.2555571 

-0.4946821 

27.5 

-3.4923413 

-3.8449325 

•3.9685697 

0.4762284 

0.1236372 

31.622777 

•3.4923413 

-3.0577627 

-3.4511702 

-0.0411711 

-0.3934075 

35.0 

-2.816404 

-3.0577627 

-3.1181619 

0.3017579 

0.0603992 

39.763536 

-2.816404 

-2.4317489 

-2.7446192 

-0.0717848 

-0.3128703 

42.5 

-2.3596907 

-2.4317489 

-2567898 

0.2082073 

0.1361491 

50.0 

-2.3596907 

-2.4317489 

-2.1827133 

-0.1769774 

-0.2490356 

Lj  norm  of  error  on  gradients : 

1.0986497 

i  0.9918611 

In  summary,  the  results  obtained  in  this  paper  confirm  that  nodal 
positioning  by  minimizing  the  stiffness  matrix  trace  leads  to  either  an 
optimum  or  near-optimum  mesh.  If  the  mesh  is  not  optimum  it  can  be  a 
starting  mesh  for  other  mesh  refinement  procedures  and  for  procedures 
using  element  division  or  element  enhancement  (h-methods  or  p-methods). 


33  AIRCRAFT  LUG  PROBLEM 

33.1  INTRODUCTION 

Accuracy  and  reliability  of  finite  element  computation  are  among  the 
most  important  considerations  in  numerical  structural  analysis.  Run  time  and 
costs  are  becoming  less  important.  Indeed,  the  costs  incurred  in  ensuring  that 
the  results  are  accurate  are  negligible  as  compared  with  the  costs  of  potential 
consequence  of  wrong  decisions  [17]. 

In  accurate  finite  element  modelling,  a  combination  of  element  size 
modification  (h-refinement)  and  element  order  modification  (p-refinement) 
provide  the  most  efficient  solution  convergence.  An  exponential  rate  of 
convergence  can  be  achieved  with  optimally  designed  meshes  and  optimal 
order  refinements,  [18,19].  From  a  practical  standpoint,  however,  it  is  often 
difficult  to  implement  and  achieve  exponential  convergence.  Nevertheless, 
Szabo  [17]  states  that  "good"  results  can  be  obtained  through  mesh  design 
along  with  element  order  refinements. 

The  most  widely  used  procedure  is  to: 
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[1]  Select  large  elements  where  the  solution  is  known  to  be  smooth. 

[2]  Select  smaller  elements  where  the  solution  is  known  to  vary  rapidly; 
as  around  points  of  singularity. 

Szabo  suggests  that  refinement  toward  the  singularity  should  be  in 
geometric  progression  with  the  ratio  of  sizes  of  about  0.1  to  0.15.  These 
are  empirical  values.  When  small  size  ratios  are  used,  bad  element  aspect 
ratios  are  often  introduced  leading  to  poor  mesh  design.  Extreme  element 
aspect  ratios  are  the  cause  for  overestimation  of  structural  stiffness  [19]. 

The  intention  of  this  study  is  to  find  a  rationale  which  will  guide 
the  analyst  in  selecting  a  good  mesh.  The  procedure  developed  herein  has 
shown  that  an  unproved  mesh  can  be  obtained  by  minimizing  the  trace  of 
the  stiffness  matrix.  The  conventional  procedure,  as  outlined  above,  is  used 
in  initial  mesh  selection.  The  mesh  is  then  improved  by  minimizing  the 
stiffness  matrix  trace  by  moving  the  grid  point  locations.  The  "improved" 
mesh  can  then  be  used  in  the  h  or  p-version  of  mesh  refinement  leading 
to  a  much  faster  convergence.  It  is  believed  that  the  effort  spent  in 
minimizing  the  trace  will  be  rewarded  in  convergence  of  the  solution  after 
fewer  h  or  p  iterations. 
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3.3.2  ANALYSIS 


Convergence  of  finite  element  solutions  are  often  measured  by  the 
value  of  the  total  strain  energy.  With  improvements  in  grid,  the  total  strain 
energy  usually  increases.  Mesh  improvements  based  on  the  total  strain 
energy  produces  convergence  in  the  average  sense.  An  analyst,  however, 
generally  wants  to  know  the  location  and  the  accurate  value  of  the  maximum 
stress  as  well  as  accurate  values  of  displacement  at  designated  points  as 
dictated  by  design  requirements. 

If  u  is  the  actual  displacement  field  and  if  u  is  the  displacement 
predicted  by  the  finite  element  method,  then  J  |u  -  u|  j  is  the  norm  of  error 
on  displacement.  Since  the  displacement  formulation  of  the  finite  element 
method  renders  the  structure  overstiff,  u  is  larger  than  u.  Therefore  u 
converges  from  below.  Hence,  from  a  practical  standpoint,  a  combined 
displacement,  stress  and  strain  energy  criterion  should  be  used. 

In  his  study  [17]  Szabo  has  used  hierarchical  basic  functions  for 
interpolation  and  has  demonstrated  p-version  convergence  in  an  analysis  of 
an  aircraft  lug  as  shown  in  Figure  3.5.  The  lug  is  0.5  inches  thick  and  the  rest 
of  the  dimensions,  as  shown,  are  in  inches.  It  is  made  of  isotropic 
material  of  modulus  of  elasticity  of  30,000  ksi  and  Poisson’s  ratio  of 
0.3.  The  lug  is  fixed  along  AB.  A  circular  pin  carrying  a  load  of  10.0  kips, 
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inclined  at  an  angle  of  45  degrees  to  the  horizontal,  fits  tightly  in  the  hole  of 
the  lug  and  exerts  pressure  on  it.  The  results  of  total  strain  energy  changes 
obtained  in  Szabo’s  modellings  are  listed  in  Tables  3.6  and  3.7  for  easy 
comparison. 

In  the  present  study,  the  standard  QUAD4,  QUAD8,  TRIA3  and 
TRIA6  elements  of  MSC  -  NASTRAN  are  used.  The  mesh  used  by  Szabo 
is  shown  in  Figure  3.6.  Note  that  elements  2  and  11  are  distorted.  In 
attempting  to  improve  the  mesh,  grid  point  8  is  moved  to  a  new  location  as 
shown  in  Figure  3.7.  Next,  grid  points  2,3,5  and  6  are  moved  to  obtain  the 
mesh  of  Figure  3.8.  Table  3.6  shows  the  effect  upon  the  total  strain  energy 
along  with  the  reduction  in  the  stiffness  matrix  trace.  (The  trace  of  the 
stiffness  matrix  is  obtained  by  using  a  series  of  DMAP  instructions  of  MSC- 
NASTRAN  as  indicated  in  the  Appendix.)  This  improvement  is  significant 
when  compared  with  results  of  Szabo  (see  Tables  3.6  and  3.7)  for  the  same 
number  of  degrees  of  freedom.  (English  units  are  employed  in  the  tables.) 
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TABLE  3.6  -  Comparison  of  Performance  of  Meshes  1,2  and  3 


Mesh  1 

Mesh  2 

Mesh  3 

Trace  (*  10°) 

7.8712 

73158 

Degrees  of  Freedom 

36 

36 

36 

Strain  Energy(*  10'2) 
Displacement  at 

1.966856 

2.068929 

2.179273 

node  12  (•  10*3) 
Von-Mises  Stress 

5.52766 

5.92351 

639045 

at  node  1 

18.15 

18.5 

1530 

at  node  4 

23.42 

24.5 

20.51 

at  node  18 

Max.Prin.  Stress 

5.663 

5.641 

5.618 

at  node  1 

-2.056 

-3.446 

-2364 

at  node  4 

25.3 

26.89 

22.24 

at  node  18 

Min.Prin.  Stress 

4.97 

4.731 

4.696 

at  node  1 

-19.09 

-19.98 

-16.34 

at  node  4 

4.361 

5.803 

4.092 

at  node  18 

-1.195 

-1-513 

-1328 

The  maximum  vertical  displacement  occurs  at  the  tip  of  the  lug  (node 
12).  In  Table  3.6  it  is  seen  that  the  lower  value  of  trace  is  associated  with  a 
higher  maximum  displacement  (node  12).  Since  finite  element  models  tend 
to  present  "stiffer  than  actual”  systems,  the  trace  minimization  demonstrates 
the  mesh  improvement.  In  mesh  2  since  only  grid  point  8  is  moved  to  lower 
the  value  of  trace,  the  stress  values  at  nodes  1  and  4  are  more  accurate 
than  those  predicted  by  mesh  1.  In  mesh  3,  elements  1  and  3  are 
larger  than  those  of  mesh  1.  Therefore,  their  centroids  are  farther 
away  from  nodes  1  and  4  which  leads  to  lower  stress  estimation  at 
those  nodes.  Thus,  an  improved  displacement  value  is  obtained  at  the 
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expense  of  stress  values.  This  indicates  that  the  mesh  near  nodes  1  and  4 
needs  to  be  refined.  Mesh  3  can  of  course  be  improved  by  further 
minimization  of  the  trace. 


TABLE  3.7  -  Variation  of  Strain  Energy  with  P*Refinement. 


Polynomial  Order 

Degrees  of  Freedom 

Strain  Energy  (*  10'2) 

1 

36 

1.73592 

2 

100 

2.76935 

3 

170 

2.83543 

4 

266 

2.85719 

5 

388 

2.86491 

6 

536 

2.86839 

7 

710 

2.86906 

8 

910 

2.86928 

Next,  mesh  1  is  refined  by  increasing  the  polynomial  order  to  2  to 
obtain  mesh  la.  Mesh  4  is  constructed  by  refining  mesh  3  using  an  h- 
version  modification  catering  to  the  regions  of  high  stress  around  nodes  1 
and  4  and  around  the  hole. 


When  singularities  are  not  present,  an  increase  in  p  will  result  in  an 
increase  in  the  rate  of  convergence.  However,  if  singularities  are  present, 
p-refinement,  with  a  given  mesh,  will  not  necessarily  result  in  an  indefinite 
increase  in  the  rate  of  convergence.  However  an  optimal  rate  of 
convergence  can  be  obtained  by  a  proper  spacing  of  the  mesh  around  the 
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singularity  [26].  Further  mesh  4  is  improved  by  reducing  the  trace  of  the 
stiffness  matrix  to  obtain  mesh  5.  Results  of  the  three  meshes  are  listed 
in  Table  3.8. 


TABLE  3.8  -  Comparison  of  Performance  of  Meshes  la,  4  and  5. 


Mesh  la 

Mesh  4 

Mesh  5 

Trace  (*  105) 

43.166 

32.577 

31.402 

Degrees  of  Freedom 

100 

96 

96 

Strain  Energy(*  10'2) 

2.994742 

2.448802 

2.440627 

Displacement  at 
node  12  (*  10'3) 

9.02065 

732149 

733578 

Von-Mises  Stress 
at  node  1 

20.14 

2334 

23.76 

at  node  4 

26.79 

30.62 

31.13 

at  node  18 

15.22 

22.45 

2035 

Max.Prin.  Stress 
at  node  1 

-5.086 

-5.737 

-5.978 

at  node  4 

29.46 

33.77 

3438 

at  node  18 

16.6 

22.27 

2034 

Min.Prin.  Stress 
at  node  1 

-22.19 

-25.68 

-26.18 

at  node  4 

6.552 

7.805 

8.1 

at  node  18 

3326 

-0344 

-0378 

Note  that  these  meshes  have  almost  the  same  number  of  degrees  of 
freedom  and  therefore  their  performances  can  be  compared  with  each 
other.  Since  quadratic  interpolation  is  is  generally  more  accurate  than 
linear  interpolation,  mesh  la  predicts  an  improved  strain  energy  value 
compared  to  the  other  two.  However,  since  the  region  around  nodes  1 
and  4  and  also  around  the  hole  have  small  elements  in  mesh  4  and  5, 
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the  stress  values  predicted  are  much  higher  than  those  of  mesh  la.  Also,  the 
displacement  values  are  improved.  The  trace  in  mesh  5  is  smaller  because 
the  elements  around  the  hole  are  less  distorted  as  compared  with  those  of 
mesh  4.  This  leads  to  a  larger  displacement  but  at  the  expense  of  a  lower 
stress  prediction  around  the  hole.  However,  since  a  smaller  trace  improves 
the  mesh  in  the  overall  sense,  the  stress  values  at  1  and  4  are  higher.  It 
should  also  be  noted  that  the  strain  energy  of  mesh  5  is  smaller  than  that  of 
mesh  4. 

The  convergence  of  finite  element  results  in  the  energy  norm  is  not 
monotonic  [28].  However,  a  larger  displacement  indicates  a  larger  work  done 
by  external  loads  and  consequently  a  lower  potential  energy.  Again  mesh  5 
could  be  improved  by  relocating  the  nodes  to  lower  the  trace. 

Next,  mesh  5a  is  constructed  by  increasing  the  polynomial  order  of 
mesh  5  to  two.  Table  3.9  can  be  used  for  a  comparative  study  of  the  three 
meshes:  mesh  1  used  by  Szabo,  mesh  la  obtained  by  a  p-extension  of  mesh 
1,  and  mesh  5a  obtained  by  combined  h  and  p-extension  along  with  the 
improvement  procedure  based  on  trace  minimization. 

Increases  in  stress  values  of  mesh  5a  from  those  of  mesh  1  vary 
from  40%  to  400%.  Also,  the  increase  in  displacement  is  65%.  This  shows 
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the  need  for  refinements  and  improvements.  A  large  strain  energy  value 
of  mesh  la  indicates  that  a  p-  version  refinement  converges  faster  on  an 
average  sense.  Higher  stress  values  of  mesh  5a  indicate  the  superiority  of 
the  h-version  refinement.  Moreover,  the  h-version  may  introduce  distorted 
elements.  To  obtain  the  best  overall  mesh  for  a  given  number  of  degrees 
of  freedom,  it  is  therefore  useful  to  improve  the  mesh  by  trace 
minimization  procedure. 

Conclusions: 


In  view  of  the  foregoing  results  a  combined  stress,  displacement  and 
strain  energy  criterion  should  be  used  to  monitor  the  convergence.  A 
combined  grid  improvement  and  refinement  procedure  should  be  used  for 
the  best  results. 

The  study  confirms  that  the  h  and  p-extensions  lead  to  improved 
meshes.  The  study  also  shows  that  the  stiffness  matrix  trace  is  a  good 
measure  of  the  quality  of  a  mesh,  especially  when  singularities  are  present. 
Therefore,  the  step  by  step  procedure  to  be  followed  by  an  analyst  is: 

[1]  Select  large  elements  where  the  solution  is  known  to  be  smooth. 
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[2]  Select  smaller  elements  where  the  solution  is  known  to  vary  rapidly. 

[3]  Improve  the  grid  by  reducing  the  stiflhess  matrix  trace. 

[4]  Perform  a  numerical  solution. 

[5]  Refine  the  grid  by  using  the  h-version  refinement  and  by 
introducing  smaller  elements  in  high  stress  areas. 

[6]  Improve  the  grid  as  in  step  3. 

[7]  Refine  the  grid  using  the  p-version  refinement. 

[8]  Continue  to  iterate  until  the  convergence  criterion  is  satisfied. 


TABLE  3.9  -  Comparison  of  Performance  of  Meshes  1,  la  and  5a. 


Mesh  1 

Mesh  la 

Mesh  5a 

Trace  (•  10s) 

8.0856 

43.166 

158.43 

Degrees  of  Freedom 

36 

100 

322 

Strain  Energy(*  10'2) 

1.966856 

2.994742 

3.020364 

Displacement  at 
node  12  (•  10*3) 

5.52766 

9.02065 

9.10035 

Von-Mises  Stress 
at  node  1 

18.15 

20.14 

25.25 

at  node  4 

23.42 

26.79 

33.71 

at  node  18 

5.663 

15.22 

23.48 

Max.Prin.  Stress 
at  node  1 

-2.056 

-5.086 

-5.897 

at  node  4 

25.3 

29.46 

36.91 

at  node  18 

4.97 

16.60 

23.90 

Min.Prin.  Stress 
at  node  1 

-19.09 

-22.19 

-27.68 

at  node  4 

4361 

6.552 

7.768 

at  node  18 

-1.195 

3.326 

-0.872 
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Figure  3.8-  Mesh  3  (Lug  Problem) 
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Figure  3.9  -  Mesh  4  (Lug  Problem) 


Figure  3.10  -  Mesh  5  (Lug  Problem) 


3.4  DISK  PROBLEM 


3.4.1  DESCRIPTION 

A  disk  with  a  uniform  thickness  of  0.5  inch  and  a  20  inch  diameter  is 
supported  at  two  points  B  and  C  on  its  perimeter.  It  is  loaded  at  a  point  A, 
on  the  perimeter  as  shown  in  Figure  3.11.  It  is  modelled  using  TRIA6 
elements  of  MSC-NASTRAN. 

As  indicated  in  the  previous  analysis  of  the  lug,  the  initial  mesh  design 
is  an  important  step  in  the  analysis.  The  circular  disk  is  axisyminetric.  If  the 
loads  are  also  axisymmetric,  then  it  is  advantageous  to  maintain  that 
symmetry  by  choosing  annular  ring  elements.  If  the  disk  is  modelled  as 
shown  in  Figure  3.12,  then  symmetry  about  four  planes  is  retained.  However, 
the  boundary  conditions  and  the  load  warrant  only  symmetry  about  one  plane 
passing  through  AD.  In  this  study,  the  interior  nodes  are  located  on  the 
circumference  of  a  circle.  Let  r„  be  the  radius  of  this  circle.  The  non- 
dimensional  parameter  £  =  rro/r0  is  varied  to  change  the  mesh  design. 

The  graph  in  Figure  3.13  shows  the  variation  of  the  trace  of  the 
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global  stiffness  matrix  and  displacement  at  the  center  and  Figure  3.14  shows 
the  strain  energy  and  displacement  under  the  load  as  a  function  of  £.  The 
trace  reaches  its  lowest  value  at  £  =  £,  =  0.53.  At  values  of  £  significantly 
different  from  £„  elements  become  distorted  leading  to  an  increase  in  the 
trace.  As  grid  points  are  moved  towards  the  point  of  application  of  the  load 
and  support,  the  modeling  of  the  region  close  to  the  periphery  improves. 
This  is  indicated  by  the  increase  in  the  values  of  displacement  under  the  load 
and  strain  energy.  However,  the  improvement  is  restricted  by  the  increasing 
distortion  of  the  elements  which  is  indicated  by  the  increase  in  trace  values. 
Therefore,  both  displacement  under  the  load  and  strain  energy  reach  their 
maximum  values  at  £  =  £,.  They  then  decrease  in  a  similar  fashion.  It  is 
evident  that  £,  and  £,  do  not  coincide.  The  point  at  the  center  of  the  disk  is 
farthest  from  the  periphery,  and  therefore,  in  accordance  with  the  theory  of 
Saint  -  Venant,  it  should  be  least  affected  by  the  load  and  the  boundary 
conditions  imposed.  The  displacement  at  the  center  reaches  the  maximum 
at  £„  the  value  at  which  the  trace  goes  to  a  minimum.  This  shows  that  a 
minimum  trace  procedure  yields  a  good  mesh  in  the  regions  away  from  the 
boundaries  and  loads.  Ho-  ever  it  is  not  the  best  mesh  for  the  specific  loads 
and  boundary  conditions  applied.  In  order  to  achieve  this,  one  has  to  use  h 
and  p  methods  of  refinements  in  the  areas  close  to  the  boundaries.  If  the 
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restriction  that  the  interior  nodes  need  to  be  on  the  circumference  of  a 
circle,  is  relaxed,  then  the  trace  minimization  procedure  would  yield  a 
mesh  that  has  the  least  element  distorsion.  Best  results  could  be  obtained 
by  iterating  on  refinement  and  improvement  steps  until  the  error  is  below 
the  tolerance  level. 

3.4.2  CONCLUSION 

This  study  shows  that  the  trace  minimization  procedure  improves  the 
mesh  in  an  overall  sense.  For  any  specific  load  and  restraint  set,  other 
mesh  refinement  techniques  need  to  be  used. 
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3.5  LAME  PROBLEM 


3.5.1  DESCRIPTION 


One  of  the  classic  problems,  used  as  a  bench  mark  by  most 
researchers,  is  that  of  a  cylinder  subjected  to  internal  or  external  pressure. 
Lame  provided  the  theoretical  solution  for  an  infinitely  long  cylinder.  Lame’s 
results  can  be  used  to  evaluate  the  accuracy  of  finite  element  results.  Since 
it  is  possible  to  have  models  of  only  finite  length,  there  is  an  inherent  error 
associated  with  the  model.  Moreover,  when  the  cylinder  is  divided  into 
elements,  discretization  errors  are  introduced.  The  focus  of  this  study  is  on 
the  minimization  of  discretization  error  by  designing  good  grid  patterns. 

Consider  a  cylinder  of  inside  radius  r0  =  5  cms,  outside  radius  r„  =  10 
cms,  and  length  L  =  40  cms  as  shown  in  Figure  3.15.  Using  symmetry  of  the 
cylinder,  only  one  half  of  the  cylinder  is  modelled  with  the  nodes  on  the  mid¬ 
section  plane  restrained  in  the  axial  direction.  Two  stacks  of 
triangular  ring  elements  of  MSC-NASTRAN  are  used  as  shown  in 
Figure  3.16.  The  nodes  common  to  both  stacks  are  arranged  to  be 
on  a  cylindrical  surface  of  radius  rB.  The  non-dimensional  parameter 
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(  =  Tm/h  is  changed  to  vary  the  mesh  pattern. 

The  graph  in  Figure  3.17  shows  the  variation  of  the  trace  of  the  global 
stiffness  matrix  and  the  variation  of  the  strain  energy  with  respect  to  £. 
Figure  3.18  shows  that  of  the  average  radial  displacement  at  the  inside 
surface.  The  average  radial  displacement  is  computed  by  adding  the  radial 
displacement  values  at  all  the  nodes  on  the  inside  surface  and  dividing  the 
sum  by  the  number  of  nodes.  The  trace  reaches  its  maximum  value 
at  £  =  -  0.7072  which  is  the  geometric  mean  of  the  inside  and  the  outside 

radii.  The  strain  energy  reaches  its  peak  at  £  =  £,.  It  is  seen  that  £,  is 
slightly  smaller  than  £,.  For  a  uniform  mesh  £  =  £u  =  0.745.  If  strain  energy 
is  used  as  a  criterion  for  convergence  of  the  Finite  element  solution,  then  the 
mesh  with  (  =  £,  would  provide  the  best  mesh.  However,  the  mesh  with 
minimum  trace  is  very  close  to  the  best  mesh.  Moreover  it  has  been  obtained 
without  solving  the  equilibrium  equations.  The  study  shows  that  the 
minimum  trace  mesh  is  an  improvement  over  the  uniform  mesh. 

Convergence  of  strain  energy  does  not  guarantee  the  convergence  of 
displacement  and  stress  values  [17].  The  average  radial  displacement  at  the 
inside  surface  reaches  its  maximum  at  £  =  £d,  where  is  smaller  than  £,. 
Once  again  the  mesh  with  minimum  trace  is  an  improvement  over  the 
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uniform  mesh  because  is  closer  to  £<*  than  £u. 


3.5.2  CONCLUSION 

The  study  confirms  that  the  trace  minimization  procedure  will 
provide  a  good  starting  mesh.  Fewer  refinement  iterations  will  be  needed 
to  achieve  convergence  in  the  finite  element  solutions. 
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Figure  3.15  *  Cylinder  (Lame  Problem) 
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Figure  3.16  -  Finite  Element  Model  of  the  Cylinder  (Lame  Problem) 
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Figure  3.18  *  Graph  of  Average  Radial  Displacement  at  the  Inside  Surface 
(Lame  Problem) 
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4.  ALGORITHM  DEVELOPMENT 


The  example  problems  have  shown  that  the  trace  minimization 
procedure  yields  either  an  optimal  mesh,  as  in  the  Prager  probiem,  or  a 
near  optimal  mesh  as  in  the  other  examples.  In  either  case  it  yields  a 
very  good  starting  mesh. 

In  any  given  problem,  it  may  not  be  difficult  to  obtain  an  expression 
for  the  trace  of  the  stiffness  matrix.  It  would,  however,  be  very  difficult  to 
obtain  recursive  relations  by  minimizing  the  trace  of  the  matrix  with 
respect  to  the  nodal  coordinates  in  order  to  obtain  the  grid  configuration. 
Instead,  let  any  arbitrary  mesh  (usually  uniform  mesh)  be  selected.  This 
mesh  may  then  be  improved  by  relocating  the  nodes  such  that  the  value 
of  the  trace  is  lowered.  The  algorithm  for  trace  minimization  is  shown  in 
the  flow  chart  in  Figure  4.1. 

There  are  three  fundamental  issues  which  are  important  for  the 
success  of  the  algorithm.  First,  the  nodes  that  should  be  relocated  need 
to  be  identified.  Second,  the  direction  and  magnitude  of  the  movement  of 
each  of  the  identified  nodes  need  to  be  determined.  Third,  a  criterion  for 
the  termination  of  the  improvement  iteration  loop  needs  to  be  established. 
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The  node  identification  step  is  relatively  simple.  Let  ku  and  k„  be 
the  largest  and  the  smallest  diagonal  entries  on  the  stiffness  matrix. 
Analogous  to  the  methods  of  bisection,  nodes  associated  with  stiffness 

entries  larger  than  1/2  (ku  +  k„)  must  be  relocated.  In  general,  the 

"cut-off'  value  could  be  expressed  as  r?c  (Ku  +  K^),  where  r;c  =  0.5  is  one 
specific  choice.  However,  for  best  results,  r)c  could  be  determined  by 
numeiical  experimentation. 

Next,  one  or  several  nodes  could  be  moved  at  a  time.  The  latter 
choice  of  the  two  will  certainly  reduce  the  CPU  time.  Let  %  be  the  set 
of  elements  connected  to  node  "i".  Let  be  the  set  of  nodes  associated 

with  elements  in  Then  the  set  xfa  can  be  described  as  the  "neighbor 

set  of  node  i".  Note  that  in  FEM,  relocation  of  node  "i"  will  effect  a 
change  in  the  stiffness  associated  with  its  neighbor  set  only.  One  of  the 
fundamental  requirements  for  better  control  in  the  process  is  to  be  able 
to  distinguish  the  effects  of  each  individual  change.  Therefore  two  nodes,  i 
and  j,  will  oualify  for  relocation  only  if  there  are  no  common  nodes  in 
their  neighbor  sets  t/>,  and  rp-y  In  mathematical  terms,  the  intersection  of 
neighbor  sets  of  all  qualifying  nodes  should  be  an  empty  set. 

The  determination  of  the  direction  and  magnitude  of  the  movement 
of  each  identified  nodes  from  its  old  location  to  its  new  location  is 
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relatively  difficult.  Observe  that  in  the  Prager  problem,  the  scalar 
coefficient  in  the  determination  of  the  trace  (Equation  3.6)  is  given  by 
EA(x)/l(x).  Any  relocation  of  the  node  which  increases  the  element  length 
also  decreases  the  cross  sectional  area.  Therefore  it  reduces  the  stiffness 
contribution  to  the  trace.  Node  relocations  can  be  accomplished  by 
observing  the  expressions  for  element  stiffnesses  and  the  role  of  each  of 
the  parameters  involved  such  as  A(x)  and  l(x)  in  the  tapered  bar  problem. 
Each  node  can  be  selected  and  moved  to  a  new  location  manually  by 
using  a  graphic  terminal.  However  the  process  is  slow,  cumbersome  and 
inefficient. 

Another  approach  is,  for  each  identified  node,  to  obtain  its  neighbor 
set.  Then  compute  the  trac~  of  the  submatrix  corresponding  to  the 
neighbor  set,  and  store  it.  The  most  important  step  in  this  approach  is 
the  determination  of  the  trace  gradient.  In  the  one  dimensional  case,  the 
gradient  can  be  computed  by  difference  formulae  once  the  value  of  the 
trace  is  known  at  another  point.  Therefore,  select  a  new  location  at  some 
distance  away.  (A  discussion  on  the  magnitude  of  this  distance  is  given  in 
the  following  paragraph).  Next,  compute  the  trace  of  the  neighbor  set 
corresponding  to  the  new  location.  The  trace  gradient  computed  will 
indicate  the  direction  and  magnitude  of  movement  for  relocation.  In  the 
two  dimensional  case,  first  compute  the  trace  of  the  submatrix 
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corresponding  to  the  neighbor  set  of  the  node  in  its  original  location  as 
described  before.  Second,  select  a  new  location  for  the  node  in  any 
direction  and  compute  the  trace.  Next,  select  another  location  in  a 
direction  perpendicular  to  the  direction  chosen  for  the  selection  of  the 
first  new  location.  Again  compute  the  trace.  Using  the  three  trace  values, 
gradients  in  the  two  mutually  perpendicular  direction  can  be  computed, 
which  when  added  vectorially  will  yield  the  gradient  at  the  original 
location.  Similarly,  in  the  three  dimensional  case,  three  new  locations  on 
three  mutually  perpendicular  lines  should  be  used.  Finally,  the  node 
should  be  moved  in  the  direction  indicated  by  the  gradient. 

The  magnitude  of  movement  for  gradient  computation  and  for  final 
node  relocation  will  be  a  fixed  percentage  of  the  distance  between  the 
node  under  consideration  and  its  neighbor  in  the  direction  of  the 
gradient.  However  the  percentage  can  be  fixed  empirically  and/or  by 
numerical  experimentation. 

The  objective  of  this  algorithm  is  to  produce  a  mesh  with  the  least 
trace  value.  Tne  hypothesis  is  that  the  trace  minimization  procedure  will 
distribute  the  stiffness  uniformly  among  the  nodes  and  elements.  Therefore 
the  criterion  for  termination  of  the  improvement  iteration  loop  can  be 
based  either  directly  on  the  decrement  in  trace  valu*’  or  the  uniformity  of 
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the  stiffness  values  at  all  node  points. 

To  see  how  uniform  stiffness  results  in  a  uniform  distribution  of 
error  and  therefore  yields  the  best  mesh,  consider  a  finite  element  model 
with  n  degrees  of  freedom  (d.o.f.).  If  the  mesh  is  to  be  refined  by 
introducing  additional  nodes,  then  it  is  necessary  to  know  the  expected 
improvement  in  error  before  a  refinement  step  is  undertaken.  O.C. 
Zienkiewicz  et.  al.  [20]  and  Peano  et.  al.  [21]  have  shown  that  if  the 
n+lth  d.o.f.  is  to  be  introduced  hierarchically,  then  the  error  in  the  energy 
norm  is: 


(4.1) 


en,l 


3  (  fin-1  "  Kn+l.nlln  ) 


K, 


n+l,n+l 


where,  fn+1  is  force  corresponding  to  the  n+lth  d.o.f.,  Kn+lin+1  is  the 
stiffness  of  the  n+lth  d.o.f.,  Kn+l  n  is  the  off-diagonal  stiffness  relating  the 
n+lth  d.o.f.  to  the  original  n  d.o.f.  system,  and  un  is  the  array  of  nodal 
displacements  of  the  n  d.o.f.  system.  The  subscripts  n,l  of  the  error  e 
refer  to  the  n  original  d.o.f.  and  the  new  d.o.f. 


Zienkiewicz  [22]  has  used  the  above  error  relation  to  define  an 
error  indicator  in  the  form: 


82 


(4.2) 
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where,  f  is  the  finite  element  residual. 


In  an  adaptive  refinement  strategy,  these  indicators  are  normally 
calculated  for  all  the  d.o.f.  corresponding  to  the  next  refinement.  The 
indicators  serve  the  purpose  of  identifying  the  region  where  refinement  is 
necessary. 


Next,  the  error  corresponding  to  the  previous  iteration  wherein  the 
nth  d.o.f.  was  added,  is: 


(43) 
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The  corresponding  error  indicator  is: 


(4.4) 
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These  derivations  are  for  the  hierarchical  finite  elements.  However, 
the  error  with  the  conventional  finite  elements  will  be  of  a  similar  form. 
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The  most  general  method  of  generating  good  grids  is  to  have  an 
equal  distribution  of  some  specified  weight  function.  (  See  Eiseman  [23] 
for  a  complete  discussion  on  adaptive  grid  generation.)  Often,  the  error  in 
the  finite  element  solution  is  used  as  the  weight  function  [24].  Therefore 
the  objective  is  to  distribute  the  error  equally  among  all  elements. 
However,  the  value  of  the  residual  f  can  be  obtained  only  after  the 

equilibrium  equations  are  solved.  Nevertheless,  one  way  of  obtaining  an 
equi-distribution  of  error  a  priori  is  by  having  uniform  element  stiffnesses. 
As  a  consequence,  f  will  be  nearly  uniform  among  the  elements.  The 

trace  minimization  procedure  developed  herein  produces  such  a  result. 

Consider  again  the  Heat  Transfer  Example  of  section  3.2.  Note  that 

each  of  the  ratios  in  the  optimality  condition,  Equation  (3.44),  can  be 

equated  to  a  constant  7. 


Substituting  into  Equation  (3.42),  the  element  stiffness  coefficient  is  : 

(4.6)  S,  =  jr  k,  i_Ll 

7  -  1 


which  is  a  constant.  Therefore  the  trace  minimization  procedure 
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produces  a  uniform  element  stiffness. 

Finally,  observe  the  graphs  of  errors  in  Figures  4.2,  4.3  and  4.4.  The 
errors  are  equally  distributed  with  the  improved  mesh.  There  is  a  skewed 
distribution  with  the  uniform  mesh.  In  order  to  compare  the  error 
distribution  among  the  elements,  rms  errors  were  calculated  using  50 
uniformly  spaced  points  along  the  length  of  each  element.  Also  the  overall 
rms  error  for  the  model  was  calculated  using  all  the  points.  Table  4.1 
shows  the  rms  error  distribution.  Note  that  in  all  the  cases,  the  improved 
mesh  distributes  the  error  more  uniformly  than  the  uniform  mesh.  The 
rms  errors  on  the  elements  are  almost  exactly  equal  in  the  case  where 
temperatures  are  specified  at  the  boundaries.  Therefore  the  mesh  obtained 
is  optimal.  Similar  results,  however,  are  not  obtained  in  the  case  where 
both  temperature  and  temperature  gradient  are  specified  because  of  the 
inability  of  FEM  to  strongly  satisfy  the  Neumann  boundary  conditions 
[25].  Nevertheless,  it  demonstrates  the  usefulness  of  the  trace  minimization 
procedure  in  a  priori  grid  refinement. 
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TABLE  4.1  -  Comparison  of  RMS  Error  Distribution  Among  the  Elements 
Between  Finite  Element  Uniform  Mesh  and  Improved  Mesh  for 
the  Two  Models  with  Different  Boundary  Conditions. 


Temperature  B.C. 

Temperature/Gradient  B.C. 

Temoerature 

Gradient 

Uniform 

Improved 

Uniform 

Improved 

Element  1 

0.5171 

0.6786 

uunn 

Element  2 

mssm 

■ 

■ 

I 

Element  3 

0.4422 

wmm 

Kmrc  jfl 

0.2392 

0.2857 

iH  1 

0.1134 

EHSH 

Overall 

0.6800 

0.5210 

f  4494 

0.3586 

0.3755 

0.3781  | 

Another  advantage  of  using  uniform  stiffness  criterion  is  that  the 
matrix  condition  number  improves.  Matrix  condition  number  x  may  be 
shown  that: 


(4.7) 

X 

where  ^  is  the  largest  and  is  the  smallest  eigen  values  of  matrix  K. 
When,  definition  (4.7)  is  used,  x  it  is  called  the  spectral  condition  number. 
If  k^  and  kn  denote  the  smallest  and  largest  diagonal  entry  of  matrix  K,  an 
expression  for  the  lower  bound  for  x  would  be: 

(4.8)  k„ 

x  >  - 

K 
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Utku  and  Melosh  [27]  have  shown  that  the  decimal  digits  lost  in 
computation  of  displacements  in  finite  elements  to  be: 

(4.9)  D  =  p  +  log,0eq  +  log19/c 

where  D  is  the  decimal  digits  lost  in  the  computation,  eq  is  the  residual 
unbalanced  error  and  p  is  the  precision  of  the  machine  in  decimal  digits. 
Note  that  for  a  mesh  with  uniform  stiffnesses,  k„  »  lq,  and  therefore  k  >  1. 
If  K  is  equal  to  unity,  then  digits  lost  in  the  computation  is  the  minimum. 
The  criterion  of  uniform  stiffness  therefore  attempts  to  obtain  results  with  the 
least  manipulation  error. 
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Figure  4.1  -  Flow  Chart  of  Trace  Minimization  Algorithm 
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Figure  4.1  (coni;) 
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Figure  4.3  -  Graph  of  Error  in  Temperature  in  the  Cylinder  Model  with  both 
Temperature  and  Temperature  Gradient  Specified  Boundary 
Conditions 


mm 


Figure  4.4  -  Graph  of  Error  in  Temperature  Gradients  in  the  Cylinder  model 
with  both  Temperature  and  Temperature  Gradient  Specified 
Boundary  Conditions 
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5.  CONCLUSIONS  AND  RECOMMENDATION 


5.1  CONCLUSIONS 

Making  a  proper  choice  of  mesh  is  an  important  step  in  the  finite 
element  analysis  for  obtaining  accurate  results.  Although  engineering 
judgment  and  prior  knowledge  of  analysis  prove  helpful  in  mesh  design, 
the  procedure  of  trace  minimization  makes  it  possible  to  obtain  good 
meshes  without  depending  upon  such  judgments  and  knowledge.  Any 
arbitrary  mesh  may  be  used.  The  procedure  then  relocates  the  nodes  such 
that  the  value  of  the  stiffness  matrix  trace  is  lowered.  As  a  consequence 
it  redistributes  the  total  error  nearly  uniformly  among  the  elements.  Thus 
the  resulting  mesh  is  either  optimal  or  near  optimal.  The  main  advantage 
of  the  procedure  is  that  a  good  mesh  can  be  obtained  before  the 
equilibrium  equations  are  solved.  A  posteriori  methods  could  be  used  to 
refine  the  mesh  even  further.  With  the  use  of  the  trace  minimization 
procedure,  fewer  a  posteriori  refinements  become  necessary  to  obtain  the 
desired  accuracy  level  than  when  the  procedure  is  not  used. 

Most  finite  element  packages  have  routines  to  check  the  correctness 
of  the  mesh  generated.  They  check  node  coincidences,  element  coincidence 
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and  element  distortions.  The  routines  iterate  on  the  elements  to  check  if 


the  distortion  is  lower  or  higher  than  a  set  tolerance.  As  demonstrated  in 
the  aircraft  lug  analysis  example,  the  trace  minimization  procedure  yields  a 
mesh  with  minimum  element  distortion  -  a  useful  byproduct.  Moreover,  if 
the  algorithm  presented  in  Chapter  4  is  implemented  properly,  the  CPU 
time  required  may  not  be  excessively  larger  than  that  required  for  the 
routines  used  for  a  distortion  check. 

Finally,  the  procedure  improves  the  stiffness  matrix  condition  number 
and  therefore  reduces  the  truncation  errors.  The  benefits  derived  obviously 
outweigh  the  cost  of  implementing  the  procedure. 


5.2  RECOMMENDATIONS 


The  algorithm  outlined  in  Chapter  4  needs  to  be  coded  for  efficient 
processing.  Some  of  the  constants,  such  as  the  cut  off  parameter,  need  to 
be  determined  by  numerical  experimentations.  The  procedure  needs  to  be 
validated  by  applying  it  with  problems  belonging  to  classes  other  than 
those  considered  in  this  dissertation  work. 

In  most  structural  problems,  the  stiffness  matrix  is  symmetric  and 
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positive  definite.  The  procedure  developed  is  based  upon  these 
assumptions.  It  will  be  useful  to  develop  similar  procedures  for  problems 

with  stiffness  matrices  that  are  non-symmetric  and  not  positive  definite. 
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APPENDIX 


MSC-NASTRAN  allows  the  user  to  compute  intermediate  values  or 
custom  build  solution  sequences  via  the  Direct  Matrix  Abstraction 
Program  known  as  DMAP  module.  These  DMAP  instructions  typically 
preceed  the  last  card  "CEND"  in  the  Executive  Control  Deck.  The 
following  set  of  DMAP  instructions  were  used  in  the  trace  calculations: 

Nastran  Executive  Control  Deck 

ID  PROBLEM3 , LUG 
TIME  1 
SOL  24 

• 

ALTER  219 

DIAGONAL  KGG/KGGD/  $ 

MATGEN  ,  /IDEN/1/NDF/0/0  $ 

DIAGONAL  IDEN/IVEC/  $ 

TRNSP  IVEC/IVECT/  $ 

SMPYAD  IVECT,KGGD, , , ,/TRACE/2////////2  $ 

MATPRN  TRACE  //  $ 

CEND 


Case  Control  Deck. 
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Where  KGG  is  the  global  stiffness  matrix,  KGGD  is  the  vector  of 
diagonal  stiffness  entries,  IDEN  is  an  identity  matrix  of  size  NDF  by 
NDF,  IVEC  is  the  vector  of  diagonal  entries  of  IDEN,  IVECT  is 
transpose  of  IVEC  and  TRACE  is  a  matrix  of  size  1  by  1  and  contains 
the  value  of  the  trace.  Note  that  in  the  instruction  for  matrix  generation 
the  variable  NDF  has  to  be  replaced  by  the  number  of  degrees  of 
freedom  in  the  model. 


102 


☆  u.S.  GOVERNMENT  PRINTING  OfTlCE:  -  74S-KI/M2I2 


NASA 

National  Aeronautics  ami 
Space  Administration 


1.  Report  No  NASA  CR-185170 

AVSCOM  TR  89-C-019 


4  Title  and  Subtitle 


Report  Documentation  Page 

I  2  Government  Accession  No. 


Mesh  Refinement  in  Finite  Element  Analysis  by  Minimization  of  the 
Stiffness  Matrix  Trace 


7  Author(s) 


3.  Recipient's  Catalog  No 


5.  Report  Date 

November  1989 


6  Performing  Organization  Code 


8  Performing  Organization  Report  No 


Madan  G.  Kittur  and  Ronald  L.  Huston 

None 

10.  Work  Unit  No. 

1L162209A47A 

505-63-51 

9.  Performing  Organization  Name  and  Address 

University  of  Cincinnati 

Department  of  Mechanical  and  Industrial  Engineering 

Cincinnati,  Ohio  45221-0072 

1 1 .  Contract  or  Grant  No. 

NSG-3188 

12.  Sponsorin'?  Agency  Name  and  Address 

Propulsion  Directorate 

U.S.  Army  Aviation  Research  and  Technology  Activity— AVSCOM 

Cleveland,  Ohio  44135-3127 
and 

NASA  Lewis  Research  Center 

Cleveland,  Ohio  44135-3191 

13  Type  of  Report  and  Period  Covered 

Contractor  Report 

Final 

14.  Sponsoring  Agency  Code 

15.  Supplementary  Notes 

Project  Manager,  Fred  B.  Oswald,  Propulsion  Systems  Division,  NASA  Lewis  Research  Center.  Madan  G.  Kittur, 
presently  employed  by  Aero  Structures,  Arlington,  Virginia.  Ronald  L.  Huston,  University  of  Cincinnati. 


16.  Abstract 

Most  Finite  Element  packages  provide  means  to  generate  meshes  automatically.  However,  the  user  is  usually  confronted  with 
the  problem  of  not  knowing  whether  the  mesh  generated  is  appropriate  for  the  problem  at  hand.  Since  the  accuracy  of  the 
Finite  Element  results  is  mesh  dependent,  mesh  selection  forms  a  very  important  step  in  the  analysis.  Indeed,  in  accurate 
analyses,  meshes  need  to  be  refined  or  re7.oned  until  the  solution  converges  to  a  value  so  that  the  error  is  below  a  predeter¬ 
mined  tolerance.  A-posteriori  methods  use  error  indicators,  developed  by  using  the  theory  of  interpolation  and  approximation 
theory,  for  mesh  refinements.  Some  use  other  criterions,  such  as  strain  energy  density  variation  and  stress  contours  for 
example,  to  obtain  near  optimal  meshes.  Although  these  methods  are  adaptive,  they  are  expensive.  Alternatively,  a-priori 
methods,  heretofore  available,  use  geometrical  parameters— for  example,  element  aspect  ratio.  Therefore,  they  are  not  adaptive 
by  nature.  In  this  study  an  adaptive  a-priori  method  is  developed.  The  criterion  is  that  the  minimization  of  the  trace  of  the 
stiffness  matrix  with  respect  to  the  nodal  coordinates,  leads  to  a  minimization  of  the  potential  energy,  and  as  a  consequence 
provide  a  good  starting  mesh  in  a  few  examples  the  method  is  shown  to  provide  the  optimal  mesh.  The  method  is  also  shown 
to  be  relatively  simple  and  amenable  to  development  ot  computer  algorithms.  When  the  procedure  is  used  in  conjunction  with 
a-posteriori  methods  of  grid  refinement,  it  is  shown  that  fewer  refinement  iterations  and  fewer  degrees  of  freedom  are  required 
for  convergence  as  opposed  to  when  the  procedure  is  not  used.  The  mesh  obtained  is  shown  to  have  uniform  distribution  of 
stiffness  among  the  nodes  and  elements  which,  as  a  consequence,  leads  to  uniform  error  distribution.  Thus  the  mesh  obtained 
meets  the  optimality  criterion  of  uniform  error  distribution. 


17  Key  Words  (Suggested  by  Author(s)) 

Finite  element 

Mesh  refinement 

Stiffness  matrix 

18.  Distribution  Statement 

Unclassified  —  Unlimited 

Subject  Category  37 

19.  Security  Classif  (of  this  report) 

Unclassified 

20  Security  Classif  (of  this  page)  21  No  of  pages  |  22.  Price* 

Unclassified  108  A06 

"For  sale  by  the  National  Technical  I-. formation  Service  Spr  n-rjMd  frvma  2.'  t6' 


NASA  FORM  IMS  OCT  86 


